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Abstract:  Design  of  efficient  Hall-effect  thrusters  (HETs)  is  strongly  dependent  on 
understanding  of  electron  transport  throughout  these  devices.  Among  the  various  mech¬ 
anisms  are  proposed  to  account  for  enhanced  electron  transport  in  the  near-field  plume 
of  HETs  are  turbulent  fluctuations  in  the  plasma.  Since  experimental  investigations  have 
shown  that  there  is  significant  energy  in  azimuthal  oscillations  of  the  plasma,  it  is  possible 
that  these  azimuthal  oscillations  in  the  presence  of  a  radial  magnetic  field  could  be  the 
cause  of  significant  axial  transport.  This  paper  presents  the  initial  development  of  a  pseu¬ 
dospectral  azimuthal-axial  hybrid-PIC  HET  code  which  is  designed  to  explicitly  resolve 
and  filter  azimuthal  fluctuations  in  the  electrostatic  potential.  Although  development  is 
still  far  from  complete,  progress  has  been  made  in  implementing  the  electrostatic  poten¬ 
tial  solver.  Through  significant  effort,  including  the  use  of  both  numerical  limiters  and  a 
smooth  fluid  neutral  description,  stable  discharges  have  been  demonstrated  and  interim 
results  are  presented. 


Nomenclature 


E  =  magnitude  of  electric  field 

B  =  magnitude  of  magnetic  field 

E  =  electric  field 

B  =  magnetic  field 

Te  =s  electron  temperature 

p  =  electron  pressure 

Tn  =  neutral  temperature 

mn  =  neutral  mass 

nn  =  neutral  density 

ne  =  electron  density 

ni  =  ion  density 

ki  =  ionization  rate 

ks  =  Boltzmann  constant 

e  =  elementary  charge 

v  =  electron  momentum  exchange  collision  frequency 
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je  =  electron  current  density 

ji  =  ion  current  density 

fi  =  scalar  electron  mobility 

/i_L  =  perpendicular  electron  mobility 

/iA  =  electron  mobility  perpendicular  to  both  E  and  B 

< j)  =  electrostatic  potential 

me  =  electron  mass 


I.  Introduction 

The  primary  challenge  in  HET  simulation,  as  with  most  magnetized  plasma  devices,  is  the  need  to  accu¬ 
rately  model  the  transport  of  electron  across  magnetic  field  lines.  The  presumed  mechanisms  for  this  type 
of  transport,  known  as  cross-field  electron  transport,  in  HETs  are  a  combination  of  volumetric  electron  mo¬ 
mentum  exchange  collisions,  electron-wall  collisions,  and  field  fluctuations.  The  first  mechanism  refers  to  the 
idea  that  when  an  electron  undergoes  a  collision,  it  can  abruptly  shift  it’s  guiding  center  and  thus  accomplish 
net  motion  against  the  magnetic  field.  These  collisions  can  include  both  electron-ion  and  electron-neutral 
collisions,  although  for  most  relevant  HET  plasma  conditions,  the  electron-neutral  collisions  dominate.  The 
second  mechanism,  electron- wall  collisions,  presumes  that  the  interaction  between  an  electron  and  the  device 
walls  is  the  equivalent  to  a  scattering  collision.  Finally,  the  third  mechanism  refers  to  the  family  of  field 
fluctuations  (both  resolved  and  unresolved)  which,  especially  in  the  presence  of  magnetic  fields  and  fluctu¬ 
ating  densities,  can  lead  to  net  transport  effects.  The  most  popular  reduced  model  in  HET  simulation  for 
unresolved  field  fluctuation  effect  is  the  Bohm  diffusion  (or  Bohm  transport)  model. 

Traditional  computational  models  for  HETs  fall  largely  into  two  categories  -  full-PIC  and  hybrid-PIC.  While 
full-PIC  is  highly  desirable  in  that  it  actually  resolves  nonequilibrium  behavior  for  the  electron  distribution, 
hybrid-PIC  is  the  far  more  popular  approach  as  it  removes  the  strict  need  to  resolve  high  frequency  plasma 
characteristics  such  as  the  electron  plasma  frequency  and  Debye  length.  A  common  feature  of  most  existing 
hybrid-PIC  codes  is  the  treatment  of  the  electrons  as  a  fluid  in  local  thermodynamic  equilibrium.  This  al¬ 
lows  for  the  formulation  of  an  electrostatic  solver,  based  on  current  conservation,  which  relies  heavily  on  the 
details  of  the  relationship  between  the  forcing  terms  (electric  and  pressure  forces)  and  the  resulting  electron 
current  -  i.e.  the  key  to  the  solver  is  the  accuracy  of  the  electron  mobility  model. 

The  most  popular  hybrid-PIC  codes  in  common  use,  HPHall,1  relies  on  this  quasi- ID  current  conserva¬ 
tion  approach  and  resolves  the  radial-axial  plane.  The  formulation  explicitly  relies  on  zero  derivatives  for 
all  quantities  in  the  azimuthal  direction.  Experimental  observations,  including  Ellison,2  McDonald3  and 
Sekerak,4  have  clearly  shown  that  at  readily  observable  physical  and  timescales  (centimeters  and  kilohertz), 
this  assumption  is  not  true.  Indeed,  estimates  by  MacDonald  indicate  that  a  significant  fraction  of  the  total 
discharge  energy  can  be  accounted  for  in  azimuthal  structures  with  fairly  small  wavenumber  (corresponding 
to  low  order  azimuthal  modes).  For  this  reason,  it  seems  promising  to  try  to  resolve  these  structures  com¬ 
putationally  to  assess  their  contribution  to  overall  electron  transport  in  HET  operation. 

This  paper  attempts  to  address  the  role  of  resolved  field  fluctuations  in  the  azimuthal  direction  as  a  means  for 
promoting  electron  transport  in  a  hybrid-HET  code.  To  study  this  problem,  an  axial- azimuthal  hybrid-PIC 
HET  code  is  under  development,  initially  with  a  reduced  set  of  physics.  In  particular,  dynamic  ionization 
is  ignored  throught  the  use  of  imposed  electron  temperature  profiles.  Pioneering  work  in  axial- azimuthal 
simulation  has  been  performed  by  Coche5  (full-PIC)  and  Lam  and  Fernandez6  (hybrid-PIC).  This  work 
follows  the  example  of  Lam  and  Fernandez  but  substitutes  a  spectral  description  in  the  azimuthal  direction 
to  more  naturally  resolve  oscillations  in  the  plasma. 

II.  Computational  Model 

The  code  developed  for  this  paper  is  a  2D  hybrid-PIC  code  which  resolves  an  axial-azimuthal  plane  along 
the  centerline  of  the  HET.  As  with  traditional  hybrid-PIC  simulations,  the  ion  description  (singly  charged 
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xenon  only)  is  with  PIC  and  the  electron  description  is  as  a  fluid  with  the  same  density  as  the  ions;  however, 
to  avoid  the  impact  of  statistical  noise  in  the  neutral  description  on  the  mobility  tensor,  the  neutrals  are 
modeled  with  a  single  temperature  (300  K)  fluid  description  rather  than  the  particle  description  commonly 
used  in  other  hybrid-PIC  HET  codes. 

A.  Electrostatic  Potential  Solver 

Crucially,  the  electrostatic  potential  solver  is  based  on  the  same  general  current  conservation  laws  as  HPHall, 
but  it  does  not  rely  on  any  corrections  to  the  electron  mobility  with  tunable  model  parameters.  The 
foundation  of  electrostatic  potential  solver  is  current  conservation  throughout  the  device  which  can  be  written 
more  as: 


V-(ji+fe)  =  0  (1) 

Assuming  that  the  time  derivative  of  the  ion  current  density  is  zero,  we  can  replace  the  divergence  of  a 
ion  current  density  with  a  scalar  ionization  term  to  get  the  overall  current  conservation  equation. 

^  ’  je  —  kionizTleTinC  (2) 

To  get  the  electron  current  density  as  a  function  of  the  electrostatic  potential  (and  other  driving  forces) , 
we  begin  with  a  Generalized  Ohm’s  Law  description  based  on  the  steady  state  electron  momentum  equation 
where  p  =  nekBTe  and  p  =  e/mev 


je  +  /L/e  x  B  =  peneE  +  pVp  (3) 

Expanding  this  results  and  performing  some  arithmentic  results  in  the  following  expression  for  the  electron 
current  density: 


je  =  0 

where: 


(4) 
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(5) 

(6) 

(7) 

(8) 


1.  Simplifications 


To  make  this  model  computationally  tractable,  a  number  of  critical  assumptions  are  now  made.  First,  it  is 
assumed  that  this  is  purely  an  axial- azimuthal  model,  so  there  is  no  conservation  equation  in  the  z-direction 
(Eq.  7).  Next,  the  magnetic  field  is  assumed  to  everywhere  exist  only  in  the  z-direction.  This  changes  Eq.  8 
from  a  3x3  to  a  2x2  matrix  with  the  following  structure: 


p  = 


1 

l+li2B2 

— liB 

1 +ijl2B2 

l^xx 

/^xy 

P- L 

~Pa 

liB 

|_i+m2b2 

1 

l+V2B2  \ 

l^yx 

tLyy_ 

P±  _ 

(9) 


Subsitution  of  the  above  components  into  Eq.  4  then  into  Eq.  2  and  some  rearrangement  leads  to  the 
following  form: 
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d  dcj).  d  dcj). 

d  dcf).  d  dcj). 

+%(',A',en'S)  +  ^(f‘±'‘en'5;)  = 

,9  dp, _ d  dp 

+  dx^^dx)  dx^^dy^ 
d  dp  d  dp 

+aiilJ^}  +  ai{lJ^ai)  +  klW 


(10) 


Applying  the  chain  rule,  moving  the  cross  terms  involved  with  the  electrostatic  potential  to  the  RHS  as 
a  source,  and  making  the  subsitutions  of  fT±  =  p±yene  and  JT^  =  //A//ene,  reveals  the  final  form: 


dln(fi_ l)  d(j) 
dx  dx 


dln(pA)d(f> 


d2(j)  d2(j) 

dx 2  dy 2 

dln(jx a)  d(\)  dln(}Ti)  dcj) 

dy  dx  dy  dy 


1  d  .  dp  1  d  .  dp  19.  dp  ,  19.  dp  , 

~  siH1 + + 

^  kinenne 


(11) 


2.  Pseudo  spectral  representation 

To  directly  capture  the  azimuthal  modes  (which  in  this  case  are  oriented  in  the  x-direction),  the  following 
approximation  is  used  for  the  electrostatic  potential  (which  also  happens  to  be  the  same  form  as  a  DFT): 


MM 


1 

N 


N-t 


k= 0 


(12) 


Where  N  *  k  =  Lx.  The  boundary  conditions  are  straightforward  since  </>(£,  0)  =  Vanode  and  (j)(x,Ly)  = 
Vcathode  implies  that  there  are  only  DC  components  (k=0  mode)  of  the  full  Fourier  expansion  at  the  anode 
and  cathode  boundaries. 


and 


Mo) 


Vanode  =  300V  if  k  =  0 
0  if  M  0 


( fik(Ly )  - 


{VCathode  —  VV 

0 


if  k  =  0 
ifk^O 


(13) 


(14) 


This  formulation  solves  for  the  Fourier  coefficients  of  the  approximation  to  <f>  rather  than  to  0  itself. 
During  an  iteration,  the  RHS  of  Eq.  11  is  evaluated  in  physical  space  then  a  FFT  is  applied  to  the  RHS  of 
Eq.  11  in  the  x-direction  at  each  z  location  to  move  the  solution  to  frequency  space.  This  provides  Ny  sets 
of  Nx  Fourier  coefficients.  The  LHS  of  Eq.  11  is  treated  as  a  finite  difference  scheme  (j  being  the  spatial 
index  in  the  y-direction)  for  the  Fourier  coefficients  which  leads  to  N  sets  of  tridiagonal  matrix  inversions  of 
the  form: 


-^2(Pk,j  + 


4*k,j- \-i  2(f)k,j  T 

Ay2 


=  RHSk,j 


(15) 
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3.  Recovering  electron  velocity 

Recovering  the  individual  velocity  components  of  the  electrons  is  fairly  simple,  since: 

je  =  —eneue  (16) 

This  leads  to  the  formulas: 


^ e,x  —  l^-L^x  l^/\Ry 


=  -fjL±tJLEx  ~ 


M-LM  dp 
ene  dx 


+  Ra  REy  + 


Ra Rdp 
ene  dy 


(17) 


and 


^ e,y  —  f^A^x  T  /a_\_Py 


=  -ii^fiEx  ~ 


Ra Rdp 
ene  dx 


-  fi±fiEy  - 


M_l Rdp 
ene  dy 


(18) 


B.  Neutral  Flow  Solver 

Neutral  advection  is  perfromed  using  a  high-order  semi-Lagrangian  advection  scheme  developed  by  Qiu.7  In 
addition  to  the  semi-Lagrangian  advection,  a  simple  Knudsen  diffusion  model,  shown  in  Eq.  19,  is  applied 
to  the  neutral  population.  In  practice,  given  the  timescale  of  the  simulation,  this  is  a  very  minor  correction 
on  the  neutral  density  and  has  minimal  impact  on  the  evolution  of  the  solution. 


dnxe 

dt 


-DKn^2nXe  ~  —VthLchannel^nXe 


(19) 


III.  Implementation 

The  basic  geometry  of  the  simulation  is  10  cm  in  the  ^-direction  (representing  the  azimuthal  channel 
dimension)  and  4.5  cm  in  the  z-direction  (representing  the  axial  channel  dimension).  The  underlying  com¬ 
putational  grid  is  divided  into  1  mm  segments  in  the  x-  and  y-directions  and  there  is  a  single  cell  of  height 
2  cm  in  the  z-direction.  The  prescribed  magnetic  field  does  not  account  for  inductive  effects  and  remains 
constant  through  the  simulation.  Further,  there  is  no  evolution  equation  for  the  electron  temperature  but 
it  maintains  the  shape  prescribed  at  the  start  of  the  simulation.  Given  this  lack  of  a  self-consistent  electron 
energy  solver,  it  is  only  necessary  to  resolve  the  ion  timescale  motion  so  a  timestep  of  50  ns  was  used  throught 
the  simulation.  The  initial  conditions  chosen  for  each  are  shown  in  Fig.  1.  The  model  used  to  evaluate  the 
ionization  rate  are  from  Appendix  E  of  Goebel.8 

1.  Particle  boundary  conditions 

A  constant  boundary  condition  of  25  mg/s  and  n  =  3el9/m3  of  Xe  at  300K  with  an  inflow  velocity  of  1911 
m/s  is  used  to  represent  neutral  injection  at  the  anode  corresponding  to  the  5  mg/s  case  described  by  Coche5 
for  a  2 cm  wide  segment  in  the  ^-direction.  Because  that  work  provides  neither  inflow  velocity  or  equivalent 
channel  thickness  for  the  2D  model,  it  is  unclear  whether  this  flow-rate  is  5x  that  of  their  work  with  the 
5mg/s  number  implied  for  the  full  annulus  rather  than  the  segment.  To  avoid  regions  of  zero  electron  density, 
a  similar  boundary  condition  is  used  to  inject  a  small  amount  of  Xe+  into  the  solution  at  anode  -  this  seed 
population  is  sampled  form  a  reservoir  with  a  density  of  1E18  1/m3,  a  temperature  of  11605  K  and  an 
inflow  drift  velocity  of  1911  m/s.  Periodic  particle  boundary  conditions  are  used  in  the  ^-direction.  In  the 
radial  direction,  specular  reflection  for  the  neutrals  and  recombination  (to  a  surface  temperature  of  300  K) 
of  the  ions  were  used.  A  novel  particle  merge/splitting  algorithm  was  also  investigated  to  constrain  the  total 
number  of  macroparticles  in  the  simulation  -  details  of  this  work  are  provided  in  a  companion  paper.9 
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Figure  1.  Imposed  electron  temperature  and  magnetic  field 


2.  Solution  methodology 

To  account  for  the  fact  that  the  electrostatic  potential  in  Eq.  11  appears  on  both  the  left-hand  side  (as 
the  solution  variable)  and  right-hand  side  (as  a  source  term)  of  the  equation,  a  fixed  point  iterative  solve 
is  used  to  converge  the  solution  at  each  iteration.  For  the  terms  in  the  right  side  of  the  equation,  simple 
finite- difference  evualutations  were  used.  In  the  future,  #- derivatives  should  also  be  calculated  in  frequency 
space  to  improve  accuracy,  but  with  the  square  spatial  cells  used  in  this  work,  improved  accuracy  in  only  a 
single  direction  would  have  a  minor  impact  on  the  overall  error  while  dramatically  increasing  the  complexity. 
Though  spectral  accuracy  is  lost,  the  pseudospectral  solve  allows  for  single-step  solution  across  space  for  a 
given  right  hand  side  enabling  rapid  convergence  of  the  fixed-point  iteration. 


Though  this  iteration  as  described  was  successful  for  non- ionizing  cases,  the  iteration  would  diverge  when 
the  dln(jjL±) / dy  term  exceeded  a  magnitude  of  approximately  1.0  as  is  the  case  in  regions  of  steep  ionization 
gradient.  To  stabilize  the  solver,  an  ad-hoc  arctangent  limiter  shown  in  Equation  was  applied  to  the  per¬ 
pendicular  log  mobility  gradient  terms  multiplying  the  electric  field  in  the  RHS.  Once  limited,  only  a  few 
iterations  were  typically  needed  to  bring  the  L2  residual  down  to  an  acceptable  tolerance  as  shown  in  the 
results  section.  The  impact  on  the  solution  of  the  limiter  will  be  assessed  in  future  work,  but  presumeably 
if  the  physical  gradients  were  sufficiently  resolved,  the  log  of  the  gradient  could  be  made  arbitrarily  small. 
Unfortunately,  the  particle  noise  is  also  inversely  related  to  the  cell  size  which  can  also  cause  the  term  to 
become  large.  The  impact  of  the  limiter  on  the  solution  is  left  to  future  work. 


dln(fjb±) 

dy 


limiter 


1 

— - — at  an 

7rA  y 


(| ln(ji±(y  +  A y)/n±(y  -  Ay ))) 


(20) 


The  code  is  single-threaded  and  requires  only  standard  desktop  computing  resources,  but  uses  the  FFTW 
compatible  interface  to  the  cuda  accelreated  cuFFT  library10  for  the  spectral  field  solve.  In  the  spectral  solve, 
both  boundary  conditions  and  the  right  hand  side  source  terms  are  de-aliased  using  the  2/3rd  dealiasing 
rule.11 


IV.  Results 

Considerable  effort  was  expended  to  demonstrate  that  this  code  could  be  run  stably  for  even  a  short 
period  of  time.  Once  the  mobility  gradient  was  limited,  the  solution  remained  stable  over  a  wide  range 
of  parameters.  Below  shows  the  fixed-point  iteration  cycle  output  for  the  300t/l  iteration  corresponding 
t  =  15/is,  as  a  representative  example  using  the  criteria  ||(<hn  —  <f>n-1||L2  <1.0e-02  shown  in  Table  1. 

Fig.  2  shows  the  time  evolution  of  the  solution.  After  an  initial  startup  transient  due  to  inconsistency  be- 
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Iteration  300: 


Iteration  301: 


Time= 

1.500005e-05  « 

NORM: 

3 . 953806e-01 

NORM: 

4.586822e-01 

NORM: 

2.932743e-01 

NORM: 

5.065295e-02 

NORM: 

1.004307e-01 

NORM: 

5.369610e-02 

NORM: 

3.419605e-02 

NORM: 

3 . 912887e-02 

NORM: 

6.203489e-03 

Time= 

1.505005e-05  « 

Table  1.  Pseudo-spectral  fixed-point  iteration  for  electrostatic  potential. 


tween  the  initial  neutral  density  and  steady  solution,  the  discharge  becomes  quite  stable  in  time.  Relatively 
little  variation  in  the  solution  after  approximately  t  =  15 fis  was  observed  in  runs  as  far  out  as  200/is.  The 
lack  of  an  observed  breathing  mode  is  attributed  to  the  constant  imposed  temperature  profile  used  in  place 
of  solving  the  electron  energy  equation,  but  confirmation  of  this  observation  is  left  to  future  work  studying 
the  full  coupled  system. 

As  the  solution  transient  clear  the  domain  after  approximately  15/is,  Fig.  3  shows  Ions  and  neutral  densities 
in  the  6  —  z  plane  where  the  anode  is  represented  by  the  bottom  surface  and  cathode  by  the  top.  The 
ionization  front  depletes  the  neutral  population,  but  never  attains  the  density  of  the  original  distribution 
because  they  are  rapidly  accelerated  out  of  the  domain  by  the  electric  field.  The  results  show  little  structure 
in  the  ^-direction  which  probably  aids  in  the  progression  towards  completely  stable  plasma  discharge.  The 
coupling  of  field  to  ion  density  alone  appears  to  not  break  the  symmetry  of  the  problem.  Once  the  energy 
equation  is  included  in  the  model,  it  will  be  interesting  to  see  if  azimuthal  structure  will  be  recovered  due 
to  the  rapid  ionization  quenching  the  electron  temperature  locally. 


Ion  Density 


Ion  Macroparticles 


Xe  Density 


1.0  2.0  3,0 

Axial  Position,  (cm) 


Number  Density 

(Mean#/m3) 


.0 

100.0  200.0  300.0  4WMJ  J.b+16 


1.0  2-0  3,0 

Axial  Position,  (cm) 


1.0  2.0  3,  D  *.0 

Axial  Position,  (cm) 


Compuational  Particles 

(Mean#  /  Cell) 


Number  Density 

(Mean#/m3) 


Figure  2.  PIC  ion  densities,  macroparticle  counts,  and  neutral  Xe  densities.  Data  is  averaged  over  0-direction  and 
plotted  over  time. 
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4.0 


Angular  Position,  (cm) 

3.000e+17 _ 3.000e+  18 _  3.000e+19 

Number  Density,  (#/m3) 

Figure  3.  Ion  (top)  and  neutral  Xe  (bottom)  densities  at  t  =  15 ys.  Mean  flow  is  from  bottom  (anode)  to  top  (cathode). 


V.  Conclusions  and  Future  Work 

These  result  represent  a  preliminary  investigation  into  a  pseudospectral  axial- azimuthal  hybrid-PIC 
methodology  for  a  Generalized  Ohm’s  law  solver.  Initial  development  required  reliance  on  numerical  artifice 
(artificial  limiters  on  some  source  terms  for  the  electrostatic  potential  solver)  to  achieve  a  stable  solution. 
Not  surprisingly,  without  a  coupled  energy  (electron  temperature)  equation,  this  effort  was  unable  to  produce 
characteristic  breathing  mode  oscillations. 

Significant  future  work  is  anticipated  for  this  effort.  In  addition  to  the  clear  need  to  incorporate  a  cou¬ 
pled  electron  energy  equation,  the  impact  the  ad-hoc  limiter  must  also  be  investigated.  It  is  unclear  whether 
this  limiter  will  continue  to  be  necessary  or  whether  a  similar  effect  can  be  accomplished  with  a  more  phys¬ 
ically  motivated  wall-loss  terms  in  the  energy  equation.  Nevertheless,  this  work  represents  an  important 
initial  step  towards  a  new  research  tool  to  better  understand  the  plasma  processes  in  HET  discharges. 
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Fig.  2.  False-color  ac -coupled  video  of  a  600-V  10-A  operating  condition  at  27  000  fps  shows  a  coherent  threefold  spoke  structure  rotating  in  the 
counterclockwise  direction..  The  frequency  of  spoke  passage  by  a  given  point  on  the  channel  is  approximately  3  kHz,  but  each  spoke's  individual  rotational 
frequency  is  only  one -third  this  value,  or  1  kHz,  corresponding  to  a  velocity  of  about  500  m/s. 
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Axial  Electron  Velocity  (m/s) 


for  new  model 


•  Based  on  experimental  observations,  it  appears  that  we 
need  only  resolve  the  first  few  (m<10)  azimuthal  modes 
to  capture  spoke  instability  -  seems  a  natural  fit  for  a 
spectral  description 

•  PIC  simulations  are  inherently  noisy  -  this  makes  it  even 
more  difficult  to  deal  with  the  stiffness  of  the  mobility 
tensor;  pseudospectral  description  would  have  the  added 
benefit  of  smoothing  the  distribution 

•  Given  time  constraints,  first  attempt  at  implementing 
pseudospectral  model  would  be  a  partial  attempt: 

—  Imposed  magnetic  field  and  electron  temperature 
profile  (actually  introduced  additional  difficulties) 
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Overall  concept 


•  Azimuthal-axial  (x-y)  computational  domain 

—  Actually  an  x-y-z  code  where  z  is  one  large  cell  (2  cm) 
across 

•  Hybrid-PIC  code 

—  Electron  are  quasineutral  fluid  characterized  by  an 
imposed  temperature  profile 

—  Ions  are  macroparticles  (PIC) 

—  Neutrals  are  a  fluid  (semi-Lagrangian  advection  and 
simple  Knudsen  diffusion  model) 

—  Electrostatic  potential  is  constructed  of  Fourier  basis 
functions  in  azimuthal  direction 
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Formulation  of  electron  mobility 


From  Generalized  Ohm’s  law: 

je  +  fife  x  B  =  aeneE  +  p57p 


Retrieve  full  Cartesian  form  of  electron  mobility 


jieneEm  //  . 


dp 


dx 


Ey  -h  /' 


=  peneEz  +  /./ 


— penc 


I 

1  +  ^2i?2  +  ^2^2  +  fl2B2 


1  +  VaB* 
ft  Bz  +  }rBxBy 

- f-t  By  +  f.l‘2  BXBZ 


-f.tBZ  +  p2BXBy  flBy  +  jj2BxBz 
1  +  //2£y  ~ftBX  +  p?  By  Bz 

f.iBr  +  f.t“  By  Bz  1  +  /./ 2  B2 
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2D  electron  mobility  (x-y) 


Assuming  an  azimuthal  axial  model  where  the  x-direction  represents 
the  azimuthal  direction  and  the  y-direction  represents  the  axial  direction 
and  assuming  that  the  B  field  is  all  in  the  z-direction 

1 
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P  xy 
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B±  _ 
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1  +  (pB)2 

Starting  from  current  conservation,  after  much  manipulation  this  is  the 
formulation  for  the  potential  solver:  d24>  02d 
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Pseudospectral  representation 


Fourier  mode  basis  functions  for  potential  in  x-direction: 


N  *  k  =  L 


X 


K=0  boundary  conditions  only  on  y-ends  of  solution: 


<M0) 


v^^mv  if  k  =  0 

0  if  k  /  0 


( L  y  J 


|  Vcathode  =  Ok  if  t  =  0 
1 0  if  k  7^  0 


Potential  solver  in  y-direction  is  simply  finite-difference  of  Fourier 
coefficients: 


<Jk.j+l  +  (pkj—l 


RHSkJ 


RHS  contains  ^so  point  nonlinear  iteration  used  to  converge  solution 
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Initial/Boundary  Conditions 


In  z-direction,  ion 
recombination  to 
form  300K  Xe 
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Periodic  in  x  Outflow 


Compromises  to  stability 


Generating  stable  solution  was  much  harder  than 
anticipated  -  numerous  numerical  tricks  were  employed  to 
achieve  stability 

-  To  avoid  regions  with  no  Xe+,  Xe+  was  injected  into  the  simulation 
from  the  anode  from  a  reservoir  with  density  1E18  1/m3,  11605  K 
and  drift  velocity  of  191 1  m/s 

-  Neutrals  were  initially  PIC  but  switched  to  fluid  description 

-  Limiter  used  to  restrain  unruly  term  on  RHS 


+  Ay)/ii±(y 


at  an 


-  Upper  1/3  of  Fourier  coefficients  are  filtered  to  prevent  aliasing 
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Azimuthal-average  profiles 
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Summary  and  Future  Work 


•  Initial  efforts  to  construct  a  hybrid-PIC  azimuthal-axial 
HET  code  with  a  pseudospectral  potential  representation 
have  met  mixed  success 

—  Achieved  stable  discharge 

—  Several  undesirable  fixes  to  stabilize  discharge 

•  Significant  future  work  is  planned 

—  Time  dependent  electron  temperature  equation 

—  Pseudospectral  representation  for  terms  evaluated  on 
the  RHS  of  the  potential  equation  /  removal  of  ad-hoc 
limiters  for  term  on  RHS 
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